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Abstract 

The response of a noisy integrate-and-fire neuron with reset to periodic input is 
investigated. We numerically obtain the first-passage-time density of the pertain- 
ing Ornstein-Uhlenbeck process and show how the power spectral density of the 
resulting spike train can be determined via Fourier transform. The neuron's output 
clearly exhibits stochastic resonance. 
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1 Introduction 



Neurons are inherently stochastic information processing devices, whence the 
study of the influence of noise on neuronal signal transmission and computa- 
tion is of great interest. Since the first evidence for the enhancement of signals 
by noise was presented about 15 years ago [1], the phenomenon of stochastic 
resonance has been demonstrated in a number of physical [2,3] and biological 
systems, especially in sensory neurons [4-6]. In the wake of these experiments, 
the theory of stochastic resonance for dynamic systems has been well devel- 
oped [7,8], and was recently extended to aperiodic signals [9]. 

As neurons in higher centers of the brain need to maintain a high signal-to- 
noise ratio as well as peripheral ones, it is plausible to presume that stochastic 
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resonance is a general principle of biological information processing. Indeed, 
models describing neurons as bistable elements have been discussed in de- 
tail [10,11]. For quantitative comparison with neurophysiological data, though, 
model neurons closer to biological reality need to be investigated. To this end, 
we study in this letter the response to a sinusoidal stimulus with superimposed 
white noise of a widely used model neuron, the leaky integrate-and-fire neuron 
which is reset upon firing [12]. In this model, the development over time of the 
membrane potential is given by the solution of the Fokker-Planck equation 
describing the overdamped limit of the Ornstein-Uhlenbeck process with an 
absorbing boundary. 

Unfortunately, no analytic solution to this boundary value problem is known [13], 
while existing approximate solutions are limited to particular parameter ranges; 
in particular, they require sufficiently strong input noise [14]. Therefore, we 
numerically solve for the first-passage-time density (FPTD), i.e. the mathe- 
matical equivalent of the inter-spike-interval distribution (ISI), using a com- 
putationally efficient integral equation approach. From the FPTD, we then 
calculate the power spectral density (PSD) of the spike train generated by the 
model neuron via fast Fourier transform, employing results from the theory 
of point processes. Finally, we determine the signal-to-noise ratio (SNR) of 
the neuron's output, which clearly exhibits stochastic resonance, i.e. SNR is 
maximal for a finite strength of input noise. 

Note the crucial difference between the model studied here and the threshold 
detector model that has been studied by several authors in recent years [15- 
17]. The former is reset after each firing, whence individual threshold crossings 
are uncorrelated and the entire spike train constitutes a renewal process. The 
latter, to which we shall refer as continuous-mode model, does not include a 
reset mechanism, but assigns one spike to each threshold crossing in positive 
direction. Thus, individual crossings are correlated and the membrane voltage 
is governed by the same Fokker-Planck equation as our model, but with natu- 
ral boundaries at ±oo, permitting analytical treatment. Indeed, Jung [18] has 
given a theory of stochastic resonance in continuous-mode threshold detectors 
based on the periodic asymptotic solution of this Fokker-Planck problem. 



2 The model 



The membrane voltage x(t) of the model neuron is governed by the Langevin 
equation of the overdamped Ornstein-Uhlenbeck process [12,19,20] 

T m x(t) = -x(t) + fi + q cos(ut + tp) + £(t) , (1) 
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where we have set the resting potential to x — 0. The membrane time-constant 
r m and the drift term p are positive constants, while q, u and tp are arbitrary 
real constants, and is Gaussian white noise with zero mean and autocor- 
relation (f (t)f (0) = 2D6 ( t - 0- The initial condition is x(0) = 0, and the 
neuron fires upon reaching the threshold voltage x(t) = x t h- both x and the 
phase of the input stimulus out + <p are reset to their values at t — 0. 

To normalize variable values, we scale time as t = t/r m and voltage as x(t) = 
x(t)/x t h, so that time constant and threshold become 1, whence p = p/x t h, 
q = q/x t h, ou = r m uj, <p — ip and D = DT m /x^ h . Thus we obtain the the 
dimensionless equation 

x(t) = -x{t) + fi + q cos(ut + £(t) , (2) 

where we have dropped the bars immediately for compactness of notation. 

As mentioned above, each approach to the threshold is independent of the 
past, because of the reset upon firing. Therefore, assuming a spike train of 
infinite duration, the firing process is a stationary renewal process [21]. We 
will solve the FPT problem in the next section before examining the spike 
train as a whole in sections 4 and 5. 



3 First-Passage-Time Density 



In this section, we present an efficient numerical method for the computation 
of the FPTD 

pit) dt = Pr {x(t) = x t h = 1 in [t, t + dt) if x(t = 0) = 0} , (3) 
the theoretical counterpart of the ISI distribution. 

The Fokker-Planck equation corresponding to the Langevin equation (2) is [22] 
d d 

—V (x,t\ x ,t ) = - -£-{-x + p + qcos(ut + p))V (x,t \ x ,t ) 

Ot OX (a \ 

Q2 l 4 J 

+ D—V{x,t\x ,t ) , 

where V (x, t \ x , t ) is the probability density that the voltage is x at time t if 
it was x at time to < t. The model is thus specified by the initial and bound- 
ary conditions V r (x, 1 1 0, 0) = S(x), V r (-oo, t \ 0, 0) = and V r (1, t \ 0, 0) = 0, 
where the index r indicates restriction to x G (—oo, 1]. No analytic solution is 
known for this boundary value problem and an approximation based on the 
method of images is valid for a limited range of parameters only [14] . 



3 



Following Schrodinger [23], we thus construct an integral equation equivalent 
to the above boundary value problem, utilizing the solution Vf (x, 1 1 Xo, to) 
of (4) for the unrestricted Ornstein-Uhlenbeck process on the entire real axis, 
i.e. with boundary conditions Vf (±00, 1 1 x ,t ) = 0. The solution is [8] 



Vf (x,t\ x ,t ) 



ft™ 2 ® 



exp 



(*-(*(*)))' 



(5) 



where the mean and variance of x(t) are (writing r\ = cot" 1 uj) 



{x{t))=n + 



VTT 



: sin(wt + ip + rj) 



_l_ e -(*-*o) 



x - p 



q 



VTT 



■. sin(cjt + V 9 + V) 



a 2 (t) — D (l — e- 2 ^*")) . 
Then, the FPTD p(t) is given by the Volterra integral equation [22] 

7^/ (1, 1 1 0, 0) = f t dsV f (l,t\l,s)p(s) . 

J 



(6) 
(7) 

(8) 



Due to the sine terms in (6), the kernel Vf (l,t \ 1, s) of the above equation 
cannot be rewritten as a function of t — s alone and a solution by Laplace 
transform is not possible. A description of the FPTD via its moments cannot 
be obtained either, since such methods are based on the Laplace transform of 
the kernel [24,25]. 

We thus solve for p(t) using standard computational techniques. Since the 
kernel has an integrable square-root singularity at t — s, we rewrite (8) as 



P/(M|0,0) = r(t)p(t) + f dsVf (1, * 1 1, s) \p(s) - p(t)} , 

Jo 



(9) 



with r(t) = Jq Vf (1,£| l,s) ds. This integral can be evaluated numerically 
and, discretizing time as tj = jh with stepsize h > 0, we obtain the following 
algorithm for calculating the FPTD [26] 



P0 = , p r , 



m = 1,2,... 



(10) 



where K mJ = Vf (l,mh \ l,jh), g m = V f (l,mh \ 0,0), r m = r(mh). p = p(0) 
follows from the initial conditions. 



The algorithm defined by (10) has proven to be stable and reliable. Over a 
wide range of parameter values, the calculated FPTDs p m are strictly non- 
negative (if numerical noise of the order of machine accuracy is excluded) 
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and the norm of the distributions approaches 1 from below as the range of 
calculation is extended towards larger t. 

We found a different integral-equation approach [27,28] to be slightly less 
stable for some interesting parameter values. In regions where both algorithms 
are stable, results agree well. 



4 Power Spectral Density 



To calculate the power spectral density (PSD) of the neuron's output, let us 
first consider a train of M 5-spikes with inter-spike-intervals Tj distributed 
according to the FPTD p{Tj): 

M m 

fM(t) = £<*(*- t m ) , t m = t j , *i = • (ii) 

m=l j=l 

Neglecting the exact shape of the spikes amounts merely to dropping a form 
factor from the spectrum, while all statistically relevant information is con- 
tained in the firing times t m (see also [18]). 

For M — > oo, this process is a stationary renewal process, the spectra of which 
have been extensively discussed in mathematical literature [29,21], where they 
are known as Bartlett spectra [30] . Applications to neuronal systems have been 
rare to our knowledge [31]. 

The one-sided power spectral density of the spike train /m(0 is given by 



S M (f2) = f M (-Q)f M (-Q) + f M (Q)f M (Q) = 2f M (Q)f M (Q) (12) 
where the bar indicates complex conjugation and 

f M (f2) = - J== f tM dtf M (t) e~ iQt (13) 

y'lTTlM JO 

is the Fourier transform. Inserting (11) and (13) into (12) yields [21] 

-i M 

S M (f2) = e-^-* fc ) 

m,k=l 

= -7- 1 + / dth M (t) e-' tm + / dt h M (t) e im 

IT t M I JO JO J 
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where we have defined 

j M-\ 

M*) = 77 ^ S(t j+k -tj-t). (15) 

j,fc=i 

Integrating /im(^) over non-overlapping intervals would give the autocorrela- 
tion histogram of the neuron firing times. In the limit of an infinite spike train, 
we obtain 

h M {t) -> h{t) and t M /M -> (r) (M — > oo) , (16) 

where /i(t) is the renewal density and (r) the mean first-passage-time. Note 
that the renewal density h(t) is not a probability density but h(t) dt is the 
probability for a spike to occur in [t, t + dt). 

From the theory of renewal processes [21] we have for Q ^ 

j- Atm ^ = _m_, (17) 

Here, p(f2) is the Fourier transform of the FPTD p(r). Performing the limit 
in (14) and inserting (17), we obtain the one-sided PSD of the infinite spike 
train 

g(A)=-Mi+ " {f2 L + («>o). as) 

1 ' 7r(r>\ l-p(^) 1 - J 1 ^ 1 ^ 

Using this result, we can compute the PSD directly from the FPTD by means 
of a discrete Fourier transform. 

For white shot noise, i.e. the Poisson process with FPTD p(r) = Aexp(— Ar), 
the terms in p(f2) in (18) cancel and a white spectrum Sp = 1/tv(t) results. 
Any deviation of S(f2) from Sp indicates the presence of a signal. For the 
Ornstein-Uhlenbeck process studied here, the spectra approach Sp quickly 
for large Q (Fig. 1). We will therefore employ Sp as the reference noise level 
in section 5. 



5 Stochastic resonance 



Having set the mathematical stage, we may now explore the response of our 
model neuron to sinusoidal input. A single parameter characterizing the input 
signal is the distance-from-threshold of the deterministic trajectory (x(t)) 

e = 1 - sup (x(t)) = 1 - ( fj, + q 
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Fig. 1. Power spectral density for p = 0.97, q = 0.03, uj = O.lir, e = 0.0014 for three 
different noise levels, corresponding to D mSLX , medium and high noise. The vertical 
dotted lines mark the input frequency to and its first harmonic. 

In denning the signal-to-noise ratio, the following difficulty arises. The reset 
mechanism introduces a second timescale into the system besides the one given 
by the input frequency. Therefore, the output spectrum instead of spikes will 
have maxima of finite width, and the locations il s of these are shifted away 
from the input frequency uo (Fig. 1). We thus search a neighborhood of the 
input frequency for the signal peak and define 



As discussed above, we use the uniform spectral density Sp of the Poissonian 
spike train with firing rate 1 / (r) as noise reference level. Note that no SNR is 
calculated if the spectrum is monotonous in [(1 — a) oo, (1 + a) uo\. 

For all data shown, we have calculated the FPTD p(t) up to t = t max such that 
Jo max P(t) d£ ^ 0-99- Unless stated otherwise below, we employed a stepsize of 
h = 0.1 and set the initial phase of the stimulus to ip — 0. Parameter sets 
for which p(t) assumed negative values were discarded unless the latter could 
clearly be identified as numerical noise. PSDs were calculated at increasing 
frequency resolutions until results became consistent. The interval width for 
searching the signal was chosen as a = 0.07. 

As the central result of our work, we show in Figs. 2 and 3 the dependence 
of the signal-to-noise ratio on the input noise strength for various values of 
drift term p, modulation amplitude q and frequency u, which correspond to 
distances-from-threshold 0.001 < e < 0.01. All data clearly show stochastic 



SNR = 



max {S(n) |(1 -a) cu <Q< (1 + a) u} 
S~p 

tt(t) max {S(f2)\ (1 - a) uo < Q < (1 + a) uo} . 



(20) 
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^=0.95, q=0.05, (O=0.2ti 

10" 6 10" 4 
D 

Fig. 2. Signal-to-noise ratio vs. input noise strength for signals with small dis- 
tance-to-threshold. From top to bottom in the legend: e = 0.0014, 0.0023, 0.0046, 
0.0077. 




D 



Fig. 3. Signal-to-noise ratio vs. input noise strength for different input frequen- 
cies but the same distance-from-threshold e = 0.0014. Here, fi = 0.95 and 
q = 0.05 x a/1 + w 2 / \fl + uS{, (p = cot" 1 uji - cot" 1 u, h = 0.1 X u>i/u, w\ = 0.1-7T. 

resonance, i.e. attain the maximal signal-to-noise ratio SNR max at a noise 
strength D m£LX > 0. 

To explain why the SNR peaks, we best turn to the properties of the deter- 
ministic solution (x(t)) and the FPTD p(t), which are shown in Fig. 4 for 
uj = 0.17T, \i = 0.97, q = 0.03 (e 0.0014); this parameter set corresponds to 
the solid line in Fig. 2 and to the spectra shown in Fig. I. For strong noise, 
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Fig. 4. Deterministic solution and first-passage-time density for \i = 0.97, q = 0.03, 
co = O.Itt and e = 0.0014 as in Fig. 1. The noise levels correspond to (from top to 
bottom): below D max , at -D m ax> at D c and at high noise. The vertical dotted lines 
mark the first and second period of the input signal. (x(t)) has been shifted for 
clarity. 

the modulation of (x(t)) becomes virtually negligible and the threshold cross- 
ing probability is concentrated in a "drift peak" at small t. This drift peak 
shifts towards t = and sharpens as the noise strength is increased, becom- 
ing similar to a T-distribution (data not shown). In the spectrum, this peak 
corresponds to a widening hump shifting towards higher frequencies (Fig. 1) 
and no signal peak is left in the vicinity of the input frequency uu. As the 
input noise strength D decreases, threshold crossings become concentrated 
around the maxima of (x(t)), and firing events are synchronized to the input 
stimulus, with the first peak of pit) dominating the distribution for D max . As 
D is reduced beyond -D maX ) the peak at the first period shrinks and the fir- 
ing probability is more evenly distributed over subsequent maxima of (x(t)). 
Therefore, a variable number of maxima is skipped before the threshold is 
reached, resulting in erratic firing and thus a decrease in SNR (Fig. 4). 

Obviously, we cannot expect stochastic resonance for e < in this system, 
for if the deterministic solution (x(t)) reaches the threshold, spikes will be 
perfectly synchronized for D — 0, although the firing frequency may be far 
from the frequency of the input signal. 

For small distances-from-threshold (e < 0.003) and low frequency (uo = 0.l7r), 
we observe stochastic resonance at very small noise strengths D mSuX , and the 
SNR decays algebraically as D is increased beyond -D max (Fig. 2). Furthermore, 
this decay exhibits a crossover between two regimes at an intermediate noise 
strength D c . For D < D C) the loss in SNR is due to the widening of the peaks 
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Fig. 5. Position of SNR maximum vs. distance- from-threshold. Data pertaining to 
identical values of [i and q but different u are connected by lines. Note that -D max 
was chosen from the set of noise strengths for which calculations were performed, 
leading to discretization effects along the ordinate. 

in the FPTD, which are located at the maxima of (x(t)), while for D > D c , 
the drift peak becomes clearly discernible, corresponding to the onset of firing 
not synchronized with the input stimulus, see Fig. 4. 



If (x(t)) remains further from thresold, either due to reduced q or increased 
uj, stochastic resonance occurs at higher input noise strengths -D max and yields 
smaller maximum values of SNR, see the lower two curves in Fig. 2. This is to 
be expected, because as the deterministic solution remains smaller, the noise 
contribution to threshold crossing must increase, reducing the synchronization 
of firing events with maxima of (x(t)). 

The input noise strength -D max at which SNR attains its maximum depends 
strongly on the distance-from-threshold e, as is demonstrated in Fig. 5. Here, 
we have plotted -D max vs. e on a double-logarithmic scale. Indeed, the location 
of the SNR maximum roughly obeys a power law -D max ~ e 7 . A least squares 
fit yields 7 ~ 1.5. The detailed dependency of -D max on e is quite complex, 
though, and not yet well understood. 

On the other hand, -D max hardly depends on the input frequency uj if the input 
amplitude q is adjusted so as to obtain the same distance-from-threshold for 
all frequencies, see Fig. 3. This behavior is to be expected from the mecha- 
nism suggested above: the maximal SNR is reached as the firing probability 
is concentrated at the maxima of (x(t)). 
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6 Conclusions and Perspectives 

In this letter, we have investigated the response of a model neuron with re- 
set mechanism to sinusoidal input with additive white noise. The inter-spike- 
interval was determined by an efficient numerical method and power spectral 
densities were obtained by exploiting the renewal properties of the spike train 
generated. These techniques permitted us to study the behavior of the model 
neuron over a wide range of parameters, especially at very low noise strengths. 
We found clear evidence for stochastic resonance, i.e. the signal-to-noise ratio 
of the neuron's output shows a distinct maximum at non- vanishing input noise. 
Further, we have proposed a mechanism underlying this effect. The results sug- 
gest that nature does indeed employ stochastic resonance to obtain optimal 
signal-to-noise ratios in an inherently noisy information processing system. A 
detailed comparison with neurophysiological data will be given elsewhere. 

In future work, two questions need to be addressed. The dependence of the 
neuron's response on the phase f of the input stimulus has yet to be studied 
in detail. We expect such work to shed more light on the detailed structure 
of the dependencies of the signal-to-noise ratio on the input noise strength 
and of the position of the SNR maximum on the distance-from-threshold. 
More importantly, though, our model shares a weakness with other studies 
of integrate-and-fire neurons [14,32]: the presumed phase reset of the input 
stimulus is not very plausible from the viewpoint of neurophysiology. Work on 
an extended model overcoming this difficulty is currently in progress. 
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